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Biological cells with all of their surface structure and complex interior stripped away are essentially 
vesicles — membranes composed of lipid bilayers which form closed sacs. Vesicles are thought to 
be relevant as models of primitive protocells, and they could have provided the ideal environment 
for pre-biotic reactions to occur. In this paper, we investigate the stochastic dynamics of a set of 
autocatalytic reactions, within a spatially bounded domain, so as to mimic a primordial cell. The 
discreteness of the constituents of the autocatalytic reactions gives rise to large sustained oscillations, 
even when the number of constituents is quite large. These oscillations are spatio-temporal in 
nature, unlike those found in previous studies, which consisted only of temporal oscillations. We 
speculate that these oscillations may have a role in seeding membrane instabilities which lead to 
vesicle division. In this way synchronization could be achieved between protocell growth and the 
reproduction rate of the constituents (the protogenetic material) in simple protocells. 

PACS numbers: 02.50.Ey, 05.40.-a, 87.16.dj 



I. INTRODUCTION 

The cell is a structural and functional unit, the build- 
ing block of any living system. Cells consist of a mem- 
brane, made of a lipid bilayer, which encloses and pro- 
tects the contents of the cell, including genetic material. 
The membrane is semi-permeable: nutrients can diffuse 
in and serve as energy to support the functioning of the 
machinery Cells undergo replication (cell division): 
this is a process by which a cell, hereafter called the par- 
ent cell, divides into two or more cells, called the daugh- 
ters. The daughter cell contains in principle an exact 
replica of the parent's inner constituents, this property 
being ultimately a prerequisite for stable living organisms 
to exist. Such a process clearly relies on the synchro- 
nization between the duplication rate of the constituents 
and the growth of the container. In modern cells this 
condition is of paramount importance and is efficiently 
realized via dedicated control mechanisms, expressed as 
pathways of nested molecular checkpoints []J . This com- 
plex and delicate machinery has evolved; presumably the 
first minimalistic cells (so-called protocells na, d 
a far more straightforward and less elaborate way of di- 
viding. So focusing on the primordial cell units postu- 
lated to be present at the origin of life on Earth, can we 
conceive of a simple, though efficient, mechanism which 
could govern the division process? A possible answer to 



this question will emerge as a result of the calculations 
carried out in this paper. 

One of the most persuasive scenarios concerning the 
origin of life on Earth identifies vesicles as protocells [ll| . 
These are tiny closed sacs in which the outer membrane 
takes the form of a lipid bilayer, and so are good candi- 
dates for a minimal cell. Despite the dramatic reduction 
in complexity as compared to modern cells, vesicles still 
display many fascina ting properties, as revealed in lab- 
oratory experiments [111 . 1 1 21 ] . They are semi-permeable 
and allow for different types of chemicals to enter the 
enclosed volume, and so sustain any reaction cycles that 
may be taking place. In addition, vesicles can grow due to 
inclusion of lipid constituents into their surface, progres- 
sively adjust their shape, and eventually divide to pro- 
duce daughter vesicles. Vesicles which are initially spher- 
ical can pass through a number of intermediate shapes 
before they divide, for instance a vesicle may first change 
into an ellipsoid, then into a dumbbell shape and finally 
into two attached spheres, at which point it will divide 
in two [121 ]. 

However it must be said there is in reality very little 
theoretical evidence that the shape of the vesicle always 
follows this particular sequence, and even less experimen- 
tal evidence. It may be more appropriate to talk about 
an ensemble of vesicles and typical pathways to the state 
where division takes place. Similarly there may only be a 
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mean time to division, although it should be noted that 
there would then be a selection process which would favor 
the types of vesicles (if they could be distinguished) which 
would undergo the division proceeding at the fastest rate. 

When modeling protocells, one needs to relate the 
mechanism of growth and division to the actual mi- 
croscopic dynamics of the internal constituents. While 
vesicles can possibly define the scaffold of prototypical 
cell models, what can one say about the internal con- 
stituents? It is customarily believed that autocatalytic 
reactions |13j might have had a role in producing complex 
molecules required for the origin of life [H|]-[lj|- A chem- 
ical reaction is called autocatalytic if one of the reaction 
products is itself a catalyst for the chemical reaction. 
Clearly the reaction will speed up as more catalyst is 
produced. If there are several catalytic reactions, rather 
than just one — an autocatalytic set [H[ — then more 
complex behavior is possible, with some reactions pro- 
ducing catalysts for other reactions. This suggests that 
the interior of the protocell might have been occupied by 
interacting families of replicators, organized in autocat- 
alytic cycles. 

Autocatalytic reactions have also been invoked in the 
context of studies on the origin of life as a possible so- 
lution of the famous Eigen paradox (l9| . This is a puz- 
zle, since it limits the size of self-replicating molecules 
to perhaps a few hundred base pairs. At odds with this 
conclusion, almost all life on Earth requires much longer 
molecules to encode their genetic information. This prob- 
lem is dealt with in living cells by the presence of enzymes 
which repair mutations, allowing the encoding molecules 
to reach sizes on the order of millions of base pairs [2(| ■ 
In primordial organisms, autocatalytic cycles might have 
provided the required degree of microscopic cooperation 
to prevent Eigen's evolutionary drive to self-destruction 
to occur. 

In this paper we will investigate the properties of au- 
tocatalytic reactions within a bounded region of space, 
which we will identify with the vesicle, the whole struc- 
ture being a reference model for a protocell. The auto- 
catalytic reactions will be taken to have the form pro- 
posed by Togashi and Kaneko [2l|, |23 |. In their work, 
Togashi and Kaneko emphasized the role played by the 
noise intrinsic to the system of elementary constituents. 
This model was recently revisited by Dauxois et al. [23| . 
who used an approach based on expanding the master 
equation in a system-size expansion [24], to make an- 
alytic progress in the description of the process. This 
approach has recently been applied to a number of pro- 
cesses in biological systems to show how large oscillations 
can emerge, sustained by the stochastic component of the 
dynamics (25|-[27]]. The analysis has also been extended 
to a spatial model (28|, and our calculations will mirror 
those in this paper. 

Therefore here we will ask what happens once space 
(i.e. microscopic particle diffusion) is incorporated into 
the model. Are the oscillations robust or, conversely, do 
they get washed out through coarse-grained averaging? 



We shall demonstrate that spatio-temporal patterns do 
emerge and influence the mass transport inside the cell. 
We will also speculate that the division of the protocell 
requires an inherent degree of synchronization which may 
be triggered by collective, spatially ordered fluctuations 
in the concentration. Building on this scenario, one can 
imagine that localized peaks in the concentration might 
develop at a given stage of the vesicle evolution. Denser 
patches could then drive an instability which could po- 
tentially lead to the distortion of the membrane and so 
to division. 



II. THE MODEL 

The model we will use is a spatial version of the auto- 
catalytic model discussed in [23| . The idea is to introduce 
a spatial coarse-graining and divide the vesicle into small 
micro-cells, within which autocatalytic reactions occur. 
The cells adjoining the membrane which forms the limit 
of the vesicle have a special status, since the membrane 
allows chemicals to diffuse in from the environment and 
out into the environment. In this paper we will focus 
only on these micro-cells — those that are adjacent to 
the boundary — and lump all the interior micro-cells 
together into an inner region. We do not give the envi- 
ronment or this inner region any spatial structure; they 
simply act as a particle reservoir for the chemicals in the 
micro-cells adjacent to the membrane. 

In each micro-cell autocatalytic reactions as specified 
in [23[ occur, see Fig. [TJ More specifically we con- 
sider k chemical species, here labeled X 3 S , with the index 
s = 1 , . . . , k labeling the species and j = 1 , f2, the f2 
micro-cells where the reactions occur. The autocatalytic 
reactions take the form [23| 

xi+xi +1 v -^2xi +1 . a) 

The reactions are taken to be cyclic, so that X J k+1 = X{. 

The spatial element of the model is introduced through 
migration of chemical species between neighboring cells. 
The boundary cells will form a periodic structure in two 
dimensions, so that a Fourier-based approach can be 
used in the analysis described below. The geometry is 
schematically depicted in Fig. [21 in a two-dimensional set- 
ting, so that the micro-cells form a one-dimensional pe- 
riodic structure. It should be emphasized that although 
the scheme is illustrated in Fig. [5] with reference to a two- 
dimensional vesicle for simplicity, the setting and analy- 
sis apply in any spatial dimension including the relevant 
three-dimensional case. If the vesicle is d+ I-dimensional, 
clearly the micro-cells will form a d-dimensional periodic 
structure. 

The migration between adjacent cells is encapsulated 
in the following relations 

X* + E ] ' X 3 ' + E 3 , (2) 
E' J + XI E j ' + XI , (3) 
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where j and f label the adjacent cells and E % represents 
the number of vacancies in cell i. We will assume that 
the capacity of each cell is N, so that sum of the number 
of molecules of each species plus the number of vacancies 
equals N for every cell. 

Finally, cell j may lose a molecule X 3 to the environ- 
ment or inner region leaving a vacancy E 3 or a gain of a 
molecule X 3 from the environment or inner region, i.e. 



To describe the model as a chemical master equation, 
we denote the number of molecules of chemical species 
s in cell j by nP s , and so the state of the system can be 
characterized by the vector n = (n 1 , n 2 , n n ) where 
n 3 = (n\, n 2 , nV). The transition rate from one state 
n', to another n, is denoted by T(n|n') — with the initial 
state being on the right. For example, the transitions 
stemming from the autocatalytic cycles are 



E> 



XI. 



(4) 



There is no need to distinguish between the environment 
and inner region; the rates 7 S and /3 S can simply be re- 
garded as the combined rates for both processes. In the 
rest of the paper we will simply refer to both these regions 
as "the environment" . 

In the following we will formulate the model in terms 
of a chemical master equation and find the mean-field 
solution as well as determining stochastic corrections to 
this which occurs when TV is finite. We will also simulate 
the stochastic dynamics and compare the results with the 
analytic formulas we obtain. 



(5) 



where within the brackets we have chosen to indicate 
only the dependence on those species which are involved 
in the reaction. The transition rates associated with the 
migration between adjacent micro-cells take the form 



T(ni-l,ni' + l\ni,ni) 
T{ni + l,n( -lK,ni') 



zQ N 



zQ N \ ^ N 

m—1 
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A ^ N 

m—1 



where z is the number of nearest neighbors that each 
micro-cell has. Finally, for the interaction with the envi- 




micro-cell j 



o x J 




FIG. 1: (Color online) The volume of the cells adjacent to 
the boundary is imagined to be partitioned into f2 micro-cells 
(see also Figure [2]). Within micro-cell j the molecular species 
interact according to the autocatalytic reactions specified by 
Eqs. JTJ). In addition, the molecules can migrate from micro- 
cell j to its nearest neighbors, e.g. micro-cell j' , as depicted 
in the cartoon. A molecule of type XI (full circle) takes over 
a vacancy (dashed empty circle) of micro-cell E 3 , and so 
transforms into X J S , leaving behind a vacancy E 3 . Finally, the 
chemical can also diffuse in from the environment, a reaction 
that in turn implies changing E 3 into X 3 . The opposite holds 
for molecules that diffuse out into the environment. 
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FIG. 2: (Color online) In the spatial autocatalytic model con- 
sidered here the vesicle is imagined to be divided into small 
micro-cells. We are specifically interested in the micro-cells 
adjoining the membrane, shown in darker outline in the fig- 
ure. These latter link up together and constitute a sort of 
inner shell, immediately adjacent to the vesicle wall. Within 
each micro-cell the chemicals interact as shown in Figure [T] 
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ronment, the transition rates are 



T(r 



11 j\ 7s K 
1K) = 7Y a : 



T(n2 + l\ni) 



CI 



K 



iV 



(7) 



In Eqs. © and (J7J, explicit use has been made of the 
condition 



^ N A ' 



(8) 



to eliminate n J E , the number of vacancies in cell j. 

The system is intrinsically stochastic and may be 
described by the probability density function, P(n,t), 
which gives the probability of finding the system in state 
n at time t. The equation which governs the dynamical 
evolution of P(n, t) is the master equation [24j], which for 
the system under consideration here takes the form 



0=1 
Q 

3 = 1 



j =1 j'ej 



(9) 



where the three terms on the right-hand side refer to the 
local terms for the chemical reactions, the migration of 
chemical species between the micro-cells, and the inter- 
action with the environment, respectively. The notation 
j' G j means that the cell j 1 is a nearest-neighbor of the 
cell j. The three terms in the master equation can be 
expressed in a concise, but transparent, form by intro- 
ducing the step operator [24| defined by 

£t,jf({<}) = f(.-.,ni±l,...), (10) 

where / is an arbitrary function. The explicit forms for 
these three terms are 



T j 

'lor 



<y33 
mig 



E ( £ ^7+iJ - 1) T (ni - 1, ni +1 + IK, ni +1 ) 

s=l 

(11) 

fc 

E (?-J\,y l) T(ni - \.„i + l\ni,n^) 



where it is understood that the operator S^j also acts 
on P(n, t) when these expressions are substituted into 
Eq. ©. In Eq. (fTTj) the cyclic nature of the reactions 
means that should be identified as n{ and Sj^i j 

should be identified as £ ^ . The explicit expressions for 
the transition rates are given by Eqs. ©-(GO). These, 
together with Eqs. (|9))- (IT3l) define the model. 

The above description is exact; no approximations 
have yet been made. At this stage we could also resort 
to direct numerical simulations of the chemical reaction 
system by use of the Gillespie algorithm [2^, [3(| . This 
method produces realizations of the stochastic dynam- 
ics which are formally equivalent to those found from 
the master equation ([5]). Averaging over many realiza- 
tions enables us to calculate quantities of interest. We 
will discuss the results of performing such simulations in 
Section IIV1 but a very accurate approximation scheme 
exists which can be used to investigate models of this 
type analytically. This is the van Kampen system-size 
expansion [24]. It is effectively an expansion in pow- 
ers of A -1 / 2 , which to leading order (A — > oo) gives 
the deterministic equations describing the system, and 
which at next-to-leading order gives finite A corrections 
to these. These latter corrections take the form of linear 
stochastic differential equations which can then be ana- 
lyzed straightforwardly, especially in the case when the 
deterministic system has approached a stable fixed point. 
The method is based on substituting the ansatz 



1 



A 



cj 



(14) 



into the master equation (JSJ . Here <j>> s (t) is the solution to 
the deterministic equation, and £jj(t) is a stochastic term 
which is the difference between the actual value vP s /N 
and 4> 3 S at time t. 

We develop this approximation in the next two sec- 
tions. In Section [IIII we carry out the analysis to leading 
order, finding the deterministic equations and the rele- 
vant fixed point. In Section HVl we carry through the cal- 
culation to next-to-leading order, investigating the linear 
stochastic differential equations by taking their Fourier 
transforms. The derivations of these equations is lengthy, 
though straightforward, and the details of the expansion 
are provided in Appendices [S] and [Bj 



III. LEADING ORDER: THE DETERMINISTIC 
EQUATIONS 



V 

' env 



E fcj'Cj - 1) W - 1, "i + IK' , nj) 



s=l 



k 



[(£ s>j -l)TK-lK) 
+ (£-}-l)T(ni + l\n^}, 



s=l 



(12) 



(13) 



In the limit where the number of molecules (includ- 
ing vacancies) in each micro-cell, A, goes to infinity, the 
system becomes deterministic and is governed by a set 
of ordinary differential equations. These are found by 
substituting the ansatz (|T4| into the master equation (|9]) 
and letting A — > oo, after the introduction of a rescaled 
time r = t/(NQ). The calculation is described in Ap- 
pendix |X1 but the same equation can also be found by 
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multiplying Eq. (J9j> by n % r and summing over all states n. 
Either way one obtains the following equation for species 
s in cell j 



dT 



+ a s (a#(1 - £ <#„) + ^ £ &<H, 

rn — 1 m—1 
k 

+ 0.(1- ^<)- 7s 0i, (15) 

m—1 

where A is the discrete Laplacian operator A/| = 
(2/z) J2j l( zj(fi — /!)• In the limit where the size of the 
micro-cells tends to zero, these equations become par- 
tial differential equations, with A becoming the familiar 
Laplacian operator. In this respect, Eq. (fTij)) general- 
izes the results of [HI to the case of a spatially-extended 
system. When turning off the migration mechanism be- 
tween neighboring micro-cells, i.e. imposing a s = for 
any species s, the spatial aspects drop out and one for- 
mally recovers the ordinary differential equations given 
in [23|. 

To proceed with the analysis, and to make contact with 
the investigation carried out in [U [23| , we shall now as- 
sign the same chemical parameters to all the species. The 
migration rate is the only exception to this, and may have 
a different value for each species. We will see later that 
this will be necessary in order to find spatio-temporal 
oscillations but also, as we will see shortly, a straightfor- 
ward analysis is still possible if we maintain in a s , and 
none of the other parameters, an explicit reference to the 
index s. We will be concerned with finding the homoge- 
neous solution of Eq. (fT5"l) . that is, the solution with no 
spatial variation. The homogeneous solution is found to 
be an attractor of the deterministic dynamics, even when 
the system is initially prepared in a non-homogeneous 
configuration. This observation follows from numerical 
simulations, but can in principle be made quantitative 
by investigating the stability of the homogeneous fixed 
point. This means that no gradient in concentration is 
allowed between neighboring micro-cells, once the asymp- 
totic regime is attained. So, when searching for fixed 
points of the dynamics, one can set the terms involving 
the Laplacian in Eq. (| 15[) to zero. Since the only depen- 
dence on s, appearing in a s , multiplies these terms, there 
is also no dependence remaining on the species type, s, 
and so the fixed points are both independent of j and of 
s. Under these conditions a unique fixed point for the 
concentration, <fi* , is easily found to be 



(16) 



The result JIBJ) is 
when dealing with the 



for any s = 1, k and j = 1, 
identical to that obtained in [23 
non-spatial homologous model. 

In [231 ], fluctuations for finite TV were shown to induce 
regular temporal oscillations in the species populations, 



so significantly altering the predicted deterministic dy- 
namics. What is going to happen in the present spatial 
context? In Section ITVl we shall investigate this, the cen- 
tral point of the paper, by focusing on the next-to-leading 
order corrections in the van Kampen expansion. 



IV. NEXT-TO-LEADING ORDER: THE 
STOCHASTIC CORRECTIONS 

Equating the terms of next-to-leading order in the 
master equation, after rescaling the time, leads to the 
Fokker-Planck equation (|B1[) which governs the proba- 
bility density function of the fluctuations. This Fokker- 
Planck equation is formally equivalent to the following 
Langevin equation [3l|, [HJ 



dr 



where 



(Xi(T)Xi(r'))=BiiS(r-T') 



(17) 



(18) 



The noise term, A|(t), in Eq. (fl7|) is Gaussian with zero 
mean and with a correlator given by Eq. (1181) , from which 
it can be seen to be white. The form of the two matrices 
M and B are discussed in Appendix [Bj They depend 
on the solution of the deterministic equation </>^(t), and 
so in principle are time-dependent, since <^ is. How- 
ever, in practice we are interested in fluctuations about 
the stationary state, 0*, and so they lose their time de- 
pendence. They also only have a non-trivial spatial de- 
pendence through the presence of the discrete Laplacian, 
because the stationary state is homogeneous. Therefore 
the calculation can be considerably simplified by taking 
the spatial Fourier transform of Eqs. (fTT)) and (|18l) . As 
discussed in Appendix IB1 this gives (see also [281 ]) 



where 



<A k (r)A k (r')> - B^a d 8^, 8{r 



(19) 



(20) 



and where k is the wavevector. Here we have as- 
sumed that the micro-cells form a hypercubic lattice in 
d— dimensions with a lattice spacing a. The matrices M k 
and B k are given by Eqs. (|B^ - (|BT3]l and Eqs. (|B15|> - 
(IB19P respectively. However the important point is that 
now k is simply a label and the matrix structure origi- 
nating from the spatial nature of the problem has been 
lost. Thus both M k and 6 k are simply k x k matrices 
(recall that k is the number of chemical species) and the 
analysis from now on is as in the non-spatial case [23| . 

As we have already stressed in this paper, fluctuations 
about the stationary state need to be taken into account, 
since they can be significant even if TV is quite large. The 
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fact we can investigate these systematically is crucially 
dependent on the linearity of Eq. p^|) and that the M k 
and S k matrices are time-independent. It means that we 
can straightforwardly take the temporal Fourier trans- 
form of Eq. (H~9l) to obtain 

k 

(-i"6 sr ~ M k ) £( w ) = A k H, (21) 

r=l 

where / denotes the temporal Fourier transform of the 
function /. Defining the matrix (— ioj5 sr — M^ r ) to be 
$ k r (w), the solution to Eq. flU) is 

k 

^M = E[* k H]^^H- ( 22 ) 

3 = 1 

From previous investigations, and the nature of the 
system, we expect that the fluctuations about the sta- 
tionary state (fll))) will oscillate, and will also be sustained 
and enhanced by a resonant effect [H, [25| . This is indeed 
what is seen. To investigate this effect systematically we 
focus our attention on the power spectrum F s (k, uS) of 
the fluctuations of species s, 

p s (k, w )^(|e k M| 2 ) = 

r— 1 u=l 

The theoretical power spectrum can be found and plot- 
ted out, for any given choice of the chemical parameters, 
from Eq. (|23|) . To make contact with earlier investiga- 
tions [23| , and aiming at elucidating the spatial effects, 
we here solely focus on the choice k = 4 and select 77 = 10, 
/? = 5/32, and 7 = 5/32. When the a s are set equal to 
zero, the communication between neighboring micro-cells 
is silenced, each spatial block behaving as an independent 
unit. Based on Eq. (|23[) , a temporal peak in the power 
spectrum is predicted to occur. The peak is approxi- 
mately located at w ~ 4, in agreement with the analy- 
sis developed in [23| . Another simple limit is when the 
a s are made equal for all of the k chemical species. In 
this case, the temporal peak gets progressively damped 
at large k, the effect being more pronounced the larger 
the values for the migration parameters. A similar phe- 
nomenon was also reported to occur in [28| . 

More interestingly, in Fig. [3l we show the theoretical 
power spectrum Eq. (|2"3")l for a s that take different values 
for each of the species in the case of a two-dimensional 
vesicle (a one-dimensional periodic lattice of micro-cells, 
i.e. d = 1). The range of variation of the a s covers sev- 
eral orders of magnitude, which in turn corresponds to 
assigning a significantly different degree of mobility to 
the species. Molecules characterized by large values of 
a s will quickly diffuse, while those with smaller a s are 
associated with relatively static, and presumably, more 
massive, species. A localized peak is clearly displayed 
suggesting that organized spatio-temporal patterns can 



spontaneously emerge, due to the inherent stochasticity 
of the system. From an inspection of Figs. 02 it is also 
evident that the power spectrum shows a clear peak for 
all four species. We found that making the a s signifi- 
cantly different among species was a simple way to pro- 
duce localized spatio-temporal patterns. We also found 
that they could be produced if (at least) one of the a s 
was sufficiently large, when compared with the others. 

The conclusion of the above analysis, as well as the ac- 
curacy of the approximations that have been employed, 
can be tested via direct numerical simulations. By aver- 
aging over many realizations, we can calculate the power 
spectra after Fourier transformation. Results of the sim- 
ulation are displayed in Fig. @] for the same choice of 
parameters as in Fig. [3J The correspondence between 
the profiles is excellent and so confirms the correctness 
of our theoretical scheme. 

In summary, we have unambiguously demonstrated 
that organized spatio-temporal cycles can emerge in a 
simple model of protocells where the constituents inside 
the vesicle interact via an autocatalytic scheme. As we 
shall argue in the following, this finding provides a pos- 
sible mechanism to drive a dynamical synchronization 
between the duplication of genetic material inside a pro- 
tocell and the division of the vesicle membrane. 



V. DISCUSSION 

In this paper we have investigated how the discrete- 
ness of the constituents in an autocatalytic chemical re- 
action can lead to spatio-temporal oscillations. The oc- 
currence of temporal oscillations in such a setting, but 
without a spatial element, has previously been stud- 
ied [21,, 23] . Similarly, such oscillations have been studied 
for a predator-prey system in a spatial framework [28|, 
but the oscillations in this case did not occur at a non- 
zero value of k. In this paper, we have combined and 
generalized these treatments, and also put them into the 
context of vesicles, which suggests an interesting conse- 
quence of the oscillatory behavior. 

We can speculate that the natural tendency of the 
chemical constituents to organize in regular spatio- 
temporal cycles can be instrumental in achieving a de- 
gree of synchronization between the outer membrane of 
the vesicle and the mixture of chemicals inside. In the 
context of protocells, these chemicals undergoing auto- 
catalytic reactions are to be interpreted as a primitive 
form of genetic material. One would expect, as a minimal 
self-consistency requirement, that within a stable popu- 
lation, a vesicle would split into two when the chemical 
material contained within it had approximately doubled 
in size. It is tempting to postulate that such a property is 
a dynamical phenomenon, the density fluctuations acting 
as a positive feedback on the vesicle growth, so signaling 
when the constituents inside the vesicle are ready for the 
splitting to take place. 

Now let us imagine that the vesicle containing the 
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FIG. 3: (Color online) Power spectra calculated from Eq. (|23[) for k — 4 species and for a two-dimensional vesicle (one- 
dimensional periodic array of micro-cells). Here N = 5000, ft = 256, r\ = 10, P = 5/32, 7 = 5/32 and a = [100, 0.001, 1, 500]. 
Each pair of panels (the three-dimensional plot and its two-dimensional projection) refers to a different chemical species. A 
localized peak is displayed predicting the existence of regular spatio-temporal patterns. 
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FIG. 4: (Color online) Numerically calculated power spectra obtained from averaging 800 realizations. Stochastic simulations 
are performed via the Gillespie algorithm. Parameters are set as in Fig. [3] 
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chemical species grows, because of the inclusion of suc- 
cessive membrane constituents from the environment in 
which it moves. Laboratory experiments indicate [llj ] 
that a vesicle filled with water or solutes is kept in a 
turgid spherical shape while growing by additional ma- 
terial of a similar kind flowing in from the outside envi- 
ronment. It is believed that the vesicle remains spherical 
until a thermodynamic instability sets in which distorts 
the structure [33|, eventually leading to fissioning. Now 
suppose that the vesicle is filled by a discrete population 
of chemical constituents, which undergo an independent 
dynamics of the autocatalytic type. As illustrated in this 
paper, the chemicals experience a first rapid evolution to- 
wards the stationary state, where enhanced oscillations 
appear due to the intrinsic finiteness of the interacting 
constituents. Such oscillations might seed an instabil- 
ity [13, HH, which could resonate with the innate abil- 
ity of the container to divide, so initiating the splitting 
process. These ideas could be extended to protocells, 
where enhanced oscillations could originate in the prim- 
itive genetic material. These oscillations could signal to 
the membrane that the genetic evolution had been virtu- 
ally taken to completion and that the fission could now 
occur, so ensuring that the genetic material is passed 
on to the daughter protocells. This is a highly specu- 
lative suggestion, which calls for further investigation in 
the context of self-consistent formulations, where both 
the membrane and the genetic material are dynamically 
evolved. 

It is clear that the work presented here can be ex- 
tended in various ways. The nature of the lattice struc- 
ture that is assumed can be generalized. For instance 
it is straightforward to include next-nearest neighbors, 
next-next nearest neighbor and so on. The analytical 
treatment is analogous, and the results the same; only 
the form of the operator Ak changes. Numerical sim- 
ulations could also be performed in higher dimensions. 
In particular, a toroidal (donut-like) cell embedded in a 
three-dimensional space can be straightforwardly simu- 
lated. The inner volume of the cell is again partitioned 
into micro-cells, and distinct diffusion rates are assigned 
to the radial and longitudinal directions. Preliminary 
simulations indicate that collective modes can develop 
giving rise to organized spatio-temporal dynamics [36} . 
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Appendix A: The van-Kampen expansion 

In this Appendix we will give more details of the appli- 
cation of the van Kampen system-size expansion to the 



master equation (|9]). A general discussion of the method 
is given in van Kampen's book [24{ and a description of 
the application to a simple model showing sustained and 
enhanced stochastic fluctuations is given in |25j . The 
calculations given below build upon those carried out for 
the no n-sp atial version of the model considered in this 
paper |23| and a spatial predator-prey model (28[. We 
will occasionally refer back to these two papers below. 

The starting point for the expansion in powers of 
N^ 1 / 2 is the ansatz flTi]). From this the following two 
results can be derived [2J|. First, the left-hand side of 
the master equation §§§ is given by 



dP(n,t) an(^,i) 



dt dt 



n k 

^'EE 

j=l s=l 



dt 



where U(C m ,t) 
£ s j may be expanded: 



1±N 



(Al) 

P{n l m ,t). Second, the step operator 

d 2 



d 



(2N)- 1 



3\2 



1± JV"'9 f i +(2N)- i d 2 J + 



(A2) 



The right-hand side of the master equation may be also 
expanded. We begin by defining new operators which are 
the coefficients of N^ 1 / 2 and A -1 in the expansion of 
the particular combinations of the step operators which 
appear in the model. These are 



1 

+ 2 



= N'^Lt 



= N~*L 



1 



+l[ N ~*{ d *i 
where the operators L± s and L?, s read 

L ls = (d d -d ii+i ), /._>- (<>, 

and where L\j and L^j read 

Ljj = i'h - d ) • L 2j = (d$ - 



2j> 



In addition 



(£s 



For each of the three terms appearing in Eq. ^ , namely 
the local, migration and environmental terms, we can 
now identify the various contributions in the van Kampen 
expansion: those of order N" 1 / 2 , those of order A^" 1 
involving a single derivative, and those of order A^ -1 but 
involving two derivatives. We will examine these in turn. 
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1. Right-hand side of the master equation: the 
TV~ 1/2 terms 



The contribution from 7^ c , defined by Eq. (fTTj). is 



Using the definition of L\ s , shifting the sum on s and 
remembering that quantities with subscripts k + 1 are to 
be identified with those with subscripts 1, we obtain 

^ = £fo.+i##+i - Vsti-itDdg , (A3) 



where the superscript (1) indicates that this is only the 
contribution to T[ J oc from terms of order TV -1 / 2 . It should 



be noted that Eq. (|A3|) operates on 

The contribution from T™ ig , defined by Eq. (fT2|) and 
using the definition of L\j, is 

s \ m m / 

To write this contribution in a way which naturally in- 
volves the Laplacian operator we add to this sum two 
terms which add to zero, namely 



o = ^E< 



>m ~ 0i£<^- 



Summing the contribution over j' € j and introducing 
the discrete Laplacian A/| = (2/z) Ylj'ej(fs ~ fs) we 



obtain 



£^C (1) = -££(a«(i-£#,)+#£a# 



(A4) 



j' 



The contribution from T<? nv , defined by Eq. dT3l) . is 
immediately found to be 



T m = Y\^_(li (j) i_P± {l _Y 



(A5) 



Adding the three terms (|A3I) - (|A5I) together, and let- 
ting them act on II(£ T l n ,£) and summing over j, we may 
equate the resulting expression to the order TV 1 / 2 term 
in Eq. (jAll) . after the rescaling of time r = t/(Nfl). 
The resulting equation describes the deterministic time- 
evolution of the species s in micro-cell j in the limit 
N — > oo, and is given by Eq. (|15[) in the main text. 



2. Right-hand side of the master equation: the TV 1 
terms with a single derivative 

These contributions are expressed in terms of the op- 
erators Lis and L\j, and so are a function of the first 



derivatives in the fluctuation variables. We proceed as 
we did for the terms of order TV" 1 / 2 . 

The contribution from T\ QC , defined by Eq. (fTTj) . is 



-7-3(2) 
'loc 



E 



Vs+1 



Li. 



= ^El^^^^-^x) 



(A6) 



where once again use has been made of the definition 
Lis, the cyclic nature of the species, and the sum on s 
has been shifted. Here the superscript (2) indicates that 
this is only the contribution to T^ oc from terms of order 
TV" 1 with a single derivative. 

The contribution from 7^f g , defined by Eq. (fT2j) . is 

^E q s p« (« (- E fm) + £ c 1 - E • 

s mm 

Performing the same manipulations as before, but also 
inserting the identities 

mm mm 

gives, after summing over j f E j, 

j' 7 s ^ s m 

- A^^a+^E A ^- A ^E^)1- ( A? ) 



Finally, the contribution from T^ nv , defined by 
Eq. (fl3|) . is immediately found to be 



(A8) 



r i(2) = V T— fi + 



3. Right-hand side of the master equation: the TV 1 
terms with two derivatives 

These terms are expressed in terms of the operator 
Li s and Lij, and so are a function of the second order 
derivatives in the fluctuation variables. We proceed as 
we did in the two previous cases. 

The contribution from Ty oc is 



7loc 3) - Q £ 7,8+1 9^2s( 



^E r ' s + i (^^+i)( 



+ : -2 



^ 



d(e s ) 2 

). (A9) 
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The contribution from 7^{ g is 

^ (3> = ^£«4[W.(i-E*£))' 



r 9 _a_ 



(A10) 



Finally, the contribution from 7^ nv is found to be 

d 2 



73g' = 5E(§<i-EA) + £ 



(All) 



Adding the six terms (|A6|) - (|A11|) together, and let- 
ting them act on Jl(^ l m ,t) and summing over j, we may 
equate the resulting expression to the order one term in 
Eq. ([AT]) , after the rescaling of time r = t/(NCl). The 
resulting equation gives the stochastic time-evolution of 
the species s in micro-cell j. It takes the form of a Fokker- 
Planck equation which we now examine. 



Appendix B: The form of the matrices M and B 

To write down the differential equation for r), it 
is convenient to combine the indices s and j, labeling the 
species and micro-cells respectively, by a single index p. 
To do this we imagine the fcfi-dimensional vector with 
components ^ as an ordered sequence of f2 vectors, each 
of k components. This can be achieved by defining p = 
(J — l)k + s, so that the component £ p represents the 
fluctuations associated with the s-th species in the j-th 
micro-cell. 

Now the terms in Eqs. (|A6I) - (|A8I) take the form of a 
single derivative of £ p acting on a linear combination of 
I = l,...,kQ and the terms in Eqs. (|A9l) - (|All|) arc 
a linear combination of second order derivatives. There- 
fore using this more compact notation the Fokker-Planck 
takes the form 



dU 



d 



i 



l.p 



d 2 u 



(Bl) 



where the matrix A can be re-written as 

I 

So to specify the Fokker-Planck equation we need to 
give the form of the two (/cfi) x (kVi) matrices M and 
B. We first note that although they do not depend on 
the fluctuations £p( r )i they do depend on the solution 
of the deterministic differential equation (fT5j) , as well as 
on the reaction rates ri s ,a s , /3 S and y s . However since we 
are only interested in the fluctuations about the station- 
ary state (f>* given by Eq, (|16[) . and since we obtained this 
solution under the assumption that r/ s , j3 s and "f s were in- 
dependent of s, we take the two matrices to only depend 



on 4>*,rj,f3,j and a s . An inspection of Eqs. (IA6|) - (IA11[) 
reveals that the only spatial dependence in M and B orig- 
inates from the discrete Laplacian. This suggests that if 
we introduce spatial Fourier transforms, we should be 
able to diagonalize M and B, and so be left with matri- 
ces only in the species space. This is most easily carried 
out by not continuing to work with the Fokker-Planck 
equation (IB1[) . but instead with the equivalent Langevin 
equation [3l|, |32[ 



dr 



where 



(\ > (T)X q (r')=B m 5(T-r'), 



(B3) 



(B4) 



and where the noise term, A p (r), in Eq. (|B3|) is Gaussian 
with zero mean. This is Eq. (I17[) in the main text, but 
using the single index notation. 

We follow the conventions and methods of [HI for the 
spatial Fourier transforms. For simplicity, we shall as- 
sume that the lattice is a d— dimensional hypercubic lat- 
tice, with lattice spacing a. Then the Fourier transform, 
/k, of a function fj, is defined by 



-ik.oj ji 



(B5) 



where we have now written the lattice site label j as a vec- 
tor to emphasize the d— dimensional nature of the trans- 
form. We may now take the spatial Fourier transform 
of the matrix M . Since the only spatial dependence 
is through the discrete Laplacian, we may decompose 
Eq. (|B2t as follows: 



A * = EE M -#' = E [ M £ s) tt + m ( j p) a£ 

y r r 

(B6) 

where the two k x k matrices M^ NS ^ and M^ SP ^ will be 
specified below. It is now straightforward to take the 
spatial Fourier transform of Eq. (IB6[) to obtain [28[ 

A^ = J2[ Mi sr S) +M^A k ]^, (B7) 



where Ak is the Fourier transform of the discrete Lapla- 
cian and is given by 



A k = - ^2 [ C0S ( k 7 a ) - 1] ■ 



(B8) 



Care should be taken not to confuse the components of 
the wavevector k, and k, the number of chemical species. 
We have denoted the 7-th component of the wave vector 
k by k 7 to help avoid this confusion. 

The quantity within the square brackets in Eq. (|BT[) 
is the spatial Fourier transform of the matrix M. It is 



12 



diagonal in k— space and so depends on the single label 
k. We may therefore write it as 



M*=M^^+M^A k , 



(B9) 



where the two matrices M^ NS ^ and M^ SP ^ may be read 
off from Eqs. (|A6[) - (|A8[) . and are given by 



(BIO) 



if r = s + 1 

ifr = s-l (Bll) 

if \s — r\> 1, 



and 



M<? p) = a. [1 + (1 - fc)^] (B12) 
= a s <f>* if s ^ r. (B13) 

The matrix M^ NP ^ is exactly the one found in the non- 
spatial version of the model |23j, which is why we have 
attached the label l NS' to it to signify the non-spatial 
contribution to M. The spatial, or 'SP', contribution is 
simply M( SP )A k . 

To take the Fourier transform of the matrix B, we note 
that out of the three terms — given by Eqs. (|A9[) - (j Al 1|) 
— from which this matrix is constructed, the only non- 
trivial spatial dependence comes from Eq. (|A10[) . We 
display the contribution containing this dependence by 
noting the following relation: 



EEf— 

Lag 



J j ej 



: EE[^ 



d - 

2 

-s. 



■■>\2 



7 J- 



<JJ'> 



(B14) 



where J<W> is equal to 1 if j' and j are nearest neighbors, 
and zero otherwise. The part of the B matrix correspond- 
ing to the expression (|B14j) is {2z5y } > — 2J < yy > ) which has 
Fourier transform 



2z — 4^ cos(k 7 a) 

7=1 



-za d A k , 



using Eq. (|B8|) and z — 2d. Therefore we may express the 
matrix B in Fourier space in a similar way to Eq. (|B9[) : 



(B15) 



where the two k x k matrices B<- NS ^ and B^ may be 
read off from Eqs. (|A9j) - ()A1 1|) . and are given by 



-a d r) ( 



/3(l-fc0*)+7</>*+2r/(0*) 2 (B16) 







if r = s + 1 
ifr = «-l (B17) 
if Is — 7*1 > 1, 



and 



/3(f p ) =-2a d a s (j)* {l-kct)*) (B18) 
Bf r p) =0 if s ^ r. (B19) 

Once again, the matrix B^ NSS> is exactly the one found in 
the non-spatial version of the model 23] , up to a factor 
of a d , which is why we have attached the label 'NS' to 
it. The spatial contribution is B^^A^. 
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